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ABSTRACT 


A technique for reversion of the cumulative binomial probability distribution and an 
efficient and rapid computer program using the technique are presented. Reversion of 
the binomial series is useful in calculating the element success probability in a specified 
parallel system with redundancy which will yield the specified system success probability 
(reliability). The program is not limited to reliability problems but can also be used in 
a variety of scientific and engineering problems involving the reversion of the binomial 
series, 
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A COMPUTER PROGRAM FOR REVERSION OF THE CUMULATIVE 
BINOMIAL PROBABILITY DISTRIBUTION 
by Dari D. Blen 
Lewis Research Center 

SUMMARY 

A technique for reversion of the cumulative binomial probability distribution is pre- 
sented. Reversion of the binomial series is useful in calculating element success proba- 
bility in a parallel system with redundancy. That is, if a given parallel system and its 
attendant success probability (reliability) are specified, calculation of the required 
element success probability is of interest. 

An efficient and rapid computer program was written in FORTRAN IV language and 
is presented along with instructions for its use. The program is not limited to reliability 
problems but can also be used in a variety of scientific and engineering problems involv- 
ing the reversion of the binomial series. 


INTRODUCTION 

The binomial probability distribution has increasing application in present day sci- 
entific and engineering problems. The cumulative binomial distribution is 


X=N 

C(X s N«|N, P) = V' — P X (1 - P) N " X 

s ' Z-s x:(n - x)i 

X=N S 

It can be used to compute success probability C of a system of N parallel elements 
where the individual element success probability P is specified along with Ng, the 
number of elements required to succeed. 

Frequently, the reversion of the cumulative binomial distribution is of interest. 
That is, if a given parallel system and its attendant success probability are specified, 



calculation of the required element success probability is desired. This type of problem 
arises, for example, in life testing procedures, sampling plans, and confidence interval 
calculations. 

In general terms, the problem is to calculate the root P of the cumulative binomial 
series, given C, N, and Ng. Various methods have been used to evaluate P. Linear 
interpolation in tables of the cumulative binomial distribution (such as refs. 1 to 6) may 
be sufficient. Under certain conditions the binomial distribution can be approximated by 
other distributions which are also tabulated. Tables of the incomplete Beta function, the 
F distribution, and the normal distribution are examples. These methods result in 
limited accuracy and/or require lengthy computation. 

A table of the reverse cumulative binomial distribution is presented in reference 7, 
but its range is limited. This report presents a computer program which rapidly cal- 
culates the root P of the general cumulative binomial series over a range of parameters 
exceeding the range of reference 7. 


SYMBOLS 


C(X > N g |N, P) 


N 

N, 

I 

P 

X 


probability of Ng or more successes, given the values of N and P 

number of elements in parallel or number of trials 

specified number of successes 

element or trial success probability 

variable number of successes 


Subscripts: 

j iteration number 

0 initial value 


EQUATIONS 

The binomial probability C of exactly Ng successes in N independent trials is 
given by 

. N , N„ N-No 

C(X = No N, P) = — P b (l - P) b (1) 

a| Ng’. (N - Ng)’. 
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where P is the probability of success for an individual trial. The cumulative binomial 
probability of Ng or more successes in N independent trials is thus 


X=N 


C(X>N S |N, P) = z: 


n: 


X=N C 


X! (N - X): 


p*(l - P) N - X 


( 2 ) 


The computer program presented herein solves for P in equation (2) for specified 
values of C, N, and Ng. 

Equation (2) specifies the upper tail of the cumulative binomial probability distribu- 
tion. The lower tail can be computed from equation (2) by the following relation: 

C(X < n s |n, P) = X - C(X £ n s |n, P) (3) 

For certain values of Ng, such as Ng < N/2, the following relations result in a 
more efficient calculation than equation (2): 

X=N 

C(X < NolN, P) = ^5 ' — (1 - P) X P N-X (4) 

b| L J X! (N - X)'. 

X=N-Ng+1 


and 


X=N 

C(X > nJn, P) = 1 - \ ' — (1 - P) X P N " X (5) 

S| L. 4 XI (N - X)’. 

X=N-N g +1 

METHOD OF SOLUTION 

The method consists of specifying C, N, and Ng in equation (2) and finding the 

value of P which satisfies the equation. An initial guess for P, designated P Q , is 

used in equation (2), and the numerical value of the summation is compared with the 

specified value of C. If the absolute value of the difference between the summation and 

£ 

the specified value of C is greater than 1x10 , a new estimate for P is used in equa- 

tion (2). This process is continued until the value of the right side of equation (2) con- 
verges to the specified left-side value. 
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The Newton-Raphson method of finding the root P was chosen because it offers 
rapid convergence and because the required derivative of the cumulative binomial series 
is simple. The (j+l) th estimate for the root of an equation from the j th estimate of 
that root with the Newton-Raphson iterative technique is stated in equation form as 


P 


i + i= p j 


F < p i> 

9F(Pj) 


( 6 ) 


(ref. 8). For the problem under consideration 


F(P) = y^_NL_pX ( 1 . P )N-X 

' XI (N - X)’. 

X=N S 


- C(X 2 : Ng|N, P) 


(7) 


and 


ffl! pV 1 (1 . P) N - N S (8) 

3P (Ng - 1)1 (N - Ng)l 

The Newton-Raphson technique requires a judicious initial guess for the root P. In 
this analysis, Pq is chosen so that convergence to the root P is guaranteed. The 
binomial series is a monotonic function; therefore, convergence will always occur if Pq 
is chosen at the point of maximum slope on the curve of C as a function of P (fig. 1). 
This point is determined by setting the second partial derivative of F(P) with respect 
to P equal to zero. When the second partial derivative 

g- F -!L) = P S (1 - P) S (No - 1) - P(N - 1)1 (9) 

ap 2 (Ng - 1)1 (N - Ng)l 

is set equal to zero, the result is 



NS ' 1 

N - 1 


( 10 ) 
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The steps in the solution for the case C(X s 6 |N = 10, P) = 0. 90 are illustrated in 
figure 1. In most cases the rapidity of the convergence by the Newton -Raphson technique 
to the desired accuracy (lxlO - ®) is such that fewer than five iterations are required. 


FORTRAN PROGRAM 

The Newton -Raphson method of solution was programmed in FORTRAN IV language 
and required a storage capacity of less than 600 36 -bit words on the Lewis IBM 7094 
computer. The three inputs to the program are C(X s NglN, P), N, and Ng; their 
FORTRAN names are, respectively, C, N, and NS (see Program Listing in the ap- 
pendix). All three inputs are placed on a single card in the following format: 

C: E10.8 
N: 14 
NS: 14 

There is no limit to the number of cards which can be used to produce output for various 
sets of input. 

The computer printout is in tabular form with five column headings: C, N, NS, P, 
and ERROR IN P. The first three columns are merely the input. The solution P is 

ii. 

printed in column four; this is the j Ln iteration on P, which satisfies the requirement 
that the series sum using Pj differ from the input C by less than 1x10"®. ERROR IN P 
is simply the absolute value of the next increment in P, that is, |Pj + i - Pj|. 

The program has been run over a wide range of C, N, and NS. For high values 
of N, intermediate steps of the calculation procedure may result in extremely small 
numerical values. If the computer used for the calculation does not signal an underflow/ 
overflow condition, erroneous results could be produced for these extreme cases. For 
values of N as large as 500, however, no difficulties were encountered. Use of the 
program for values of input outside this range, therefore, should be made with an aware- 
ness of the underflow/overflow problem. 

The program sums the cumulative binomial series from the mean term (largest 
term) outward in each direction until each additive term drops below 10" . This pro- 

cedure saves the time required to compute those terms in the tail of the distribution 
which have a negligible effect on the cumulative sum and avoids the underflow problem 
encountered when computing a term like P N where P may equal 0. 5 and N may equal 
500. 
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The FORTRAN IV listing and a sample set of output are given in the appendix. The 
running time for this sample, including input/output time, was 0. 10 minute on the Lewis 
IBM 7094 computer. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, April 17, 1968, 

120-33-05-02-22. 
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APPENDIX - FORTRAN PROGRAM LISTING AND SAMPLE OUTPUT 

Program Listing 


JANUARY 24, 1968 

REVBIN 

THIS PROGRAM USES THE NEWTON— RAPHSCN TECHNIQUE TO REVERSE THE BINOMIAL 

SERIES. INPUTS ARE C,N,NS AND THE OUTPUT IS P IN THE FOLLOWING EQUATION! 

X=N N X N-X 

C = SUM C P (1-P) 

X=NS X 

PROGRAM CEVELOPEC BY DARL BIEN 1/24/68 AT NASA-LEWIS RESEARCH CENTER 

WRITE (6,31) 

31 FORMAT ( 1H1,49X,21HX±N N X N-X /4fcX,22HC = SUM C P (1 

1-P) / 50X , 9HX=NS X // , 35X , 1 HC , 14X,1HN,9X»2HNS»12X,1HP,12X, 

210HERRUR IN P//) 

1=0 
M = 0 
K = 0 

1 READ (5*101) C , N , N S 
101 FORMAT <E10.8,2I4) 

AN = N 
ANS = NS 
IF ( N— NS 11,50,51 

50 P=C**(1./AN> 

DELP =. 0. 

GO TC 40 

51 IF (NS-1) 1,52, 1C 

52 P=1.-(1.-C)**(1./AN) 

HELP = 0. 

GO TO 40 

C INITIAL GUESS ON P SUCH THAT SECOND DERIVATIVE OF C WRT P EQUALS 0 

10 P = IANS - 1. )/( AN - 1. ) 

11 Q =. 1. - P 
NLIMIT = N - NS 
ANLIMT = NLIMIT 
NL 1=N— 1 

C 

2 NP = AN*P 
ANP = NP 
NLNP = N-NP 
ANLNP = NLNP 
NSLNP=NS— NP 
NPLNS=NP— NS 
TERM1 = 1.0 

C CALCULATION OF MEAN (NP) TERM 

DO 6 K1 = 1 , NLNP 
AK 1= K 1 

6 TERM 1 = TERM 1*( ANP +AK.1 )/A'K 1*P«* ( ANP/ ANLNP ) *Q 

C IF NP LESS THAN NS, START SERIES WITH NS TERM INSTEAD 

IF(NP-NS) 15,9,9 

15 DO 17 K 1 1 = 1 , NSLNP 
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AK 1 1 = K 1 1 

17 TERMl=TERMl*t ANLNP-AK 11 + 1* ) / ( ANP + AK 1 1 ) *P /Q 
9 CGUESS = T ERM 1 
TERM2 = TERM 1 
DR CP = 1.0 F— 12 
C 

IF { NP-NS ) 29,29,43 

C NS+1 TERM, NS +2 TERM, N TERM 

29 DO 37 K21 * NS,NL1 
AK2L =. K21 

TERM2=TERM2* ( AN- AK2 1 ) / ( AK2 1+ !• 0 ) *P /Q 
IF { TERM 2-OROP ) 1A,1A,37 

37 CGCESS =CGU ES S + TERM2 
1A GO TC 16 

C TERMS GREATER THAN NP 

4? DO 7 K2 = N P , N L l 
AK 2 = K2 

TERK2=TERM2* ( AN-AK2 > / ( AK2+ 1. 0) *P/Q 
IF { TERM 2-DRO P ) 22,22,7 

7 CGUESS = CGUESS + TERM2 

C TERMS LESS THAN NP (INCLUDING NS TERM) 

2? TERM 3 = T ERM 1 

DO 8 K 3= 1 1 N°L NS 
AK 3 = K 3 

TERM1=TERM3*< ANP-AK3+1. ) / ( ANLNP + AK3 ) *Q/P 
IF (TERM 3- DR UP) 16,16,8 

8 CGUESS = CGUESS + TERM3 

C CALCULATION OF DERIVATIVE OF BINOMIAL SERIES 

16 DERIV = /'NS 

DO A J= 1 , NL I M I T 
A J = J 

A DER I V-DER I V*( ANS+ A. ) / AJ*P** < (ANS— 1.) /ANLIMT)*Q 
THE FORM OF NEWTGN-R APH SON SOLUTION IS P2 = P i-Y ( P 1 )/ Y M PI ) 

DELP IS THE CUR RECTI UN TO BE MADE ON P IF ACCURACY OF C IS NOT SATISFIED 

DELP = CGUESS/OE ' I V - C/DERIV 
ACCRYC = ABS ( CGU SS/C-L.) 

RELATIVE ERROR 0) THE FORM (CGUESS-O/C LS COMPUTED ABOVE AS ACCRYC. IF 

LESS THAN 1.0 E-6, PRINT RESULT — IF NOT GET NEW GUESS FOR P AND REPEAT 

IF ( ACCRYC— *000 OC 1 ) AO » 5 , 6 
*5 P o p - CELP 
Q = 1. - P 
GO TC 2 

A r ERRC-K = ABS( DELP) 

I = I + i 
Ml = M+L 

IF (1-5) 206,205,201 
201 WRITE (6,21) 

21 FORMAT (/) 

M = Ml 
I =i 1 

205 IF (K) 210*210,203 
21C IF (M- 53)200,200, 110 
203 IT ( M— 57 )200, 200, 1 10 
11D WRITE i 6* A1 ) 

A 1 FORMAT! IF 1 , 3AX, 1FC, 1AX, 1HN, 9X , 2HN S * 1 2 X , 1HP , 1 2 X , 10H ERROR IN P// ) 

I = 1 
M = 0 
K = K*1 

200 WRITE (6,20) C, N , NS , P , ERROR 
20 FORMAT ( 30X , F 1 0 . 8 , 7X * 1 4 ,7X , l A , 8X , F8.6 ,8X , E 10.2 / ) 

C READ ANUTFER S FT CF DATA 

GO TC l 
FIND 
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Sample Output 



X-N N 

X 

N-X 



C = SUM C 

P (1- 

P> 



X~NS X 




c 

N 

NS 

P 

ERROR IN P 

0.9CCC0000 

10 

6 

0.732682 

0.30E-07 

0 4 90OCCOOO 

10 

7 

0.812438 

0.75E-08 

0 *90000000 

10 

8 

0.884175 

0.37E-08 

0.90000000 

10 

9 

0 .945471 

C.11E-07 

0.9CCCC00O 

10 

10 

0.989519 

0 . 

0 .90000000 

50 

30 

0.675475 

0.26E-07 

0. 90000000 

50 

35 

0.767509 

0.19E-07 

0 4 90000000 

50 

40 

0.855019 

0.56E-08 

0.90000000 

50 

45 

0.935740 

0 . 

0.90000000 

50 

50 

0.957895 

0 . 

0 • 90000000 

ICO 

60 

0.656162 

0.19E-08 

0.90000000 

ICO 

70 

0.750898 

0.19E-07 

0 *90000000 

ICO 

80 

0.842332 

0 • 1 IE— 07 

0 .90000000 

ICO 

90 

0.928702 

0.10E-07 

0.90000000 

ICO 

100 

0*998947 

0 . 

0.90000000 

5 CO 

300 

0.626771 

0.14E-07 

0.90000000 

5C0 

350 

0.724691 

0.67E-07 

0.90000000 

5C0 

400 

0.821096 

0.85E-07 

0.90000000 

5C0 

450 

0.915134 

0.18E-07 

0.90000000 

5C0 

500 

0.959789 

0 . 

0.93000000 

10 

6 

0.777559 

0.60E-07 

0.95000000 

10 

7 

0.849971 

0.55E-06 

0.95000000 

10 

8 

0.912736 

0.37E-07 

0 • 9 5C CO 000 

10 

9 

0.963229 

0 . 

0 i 950 00000 

10 

10 

0.994884 

0 . 

0.95000000 

50 

30 

0.658616 

0 .67 E— 07 

0.95000000 

50 

35 

0.787896 

0.60E-07 

0 • 95000000 

50 

40 

0.871443 

0.37E-08 

0 i 95000000 

50 

45 

0.946429 

0.10E-06 

0 • 95000000 

50 

50 

0.998975 

0 . 

0.95000000 

ICO 

60 

0.673014 

0 # 16E— 06 

0.95000000 

ICO 

70 

0.765991 

0.60E-07 

0 . 95000000 

ICO 

80 

0.854755 

0.26E-07 

0.95000000 

100 

90 

0.937075 

0.19E-08 

0.95000000 

100 

100 

0.959487 

0 . 

0 J 9 5000000 

5C0 

300 

0.634565 

0.56E-08 

0.95000000 

5C0 

350 

0.731838 

0.17E-06 

0.95000000 

5C0 

400 

0.827167 

0 • 63E— 07 

0.95000000 

5C0 

450 

0.919465 

0.18E-07 

0.95000000 

5C0 

500 

0.959897 

0 . 
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Figure 1. - Element success probability and cumulative probability for six or more 
successes among ten elements. 
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